Fluctuating surface-current formulation of radiative heat transfer: theory and applications 
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We describe a novel fluctuating- surface current formulation of radiative heat transfer between bodies of ar- 
bitrary shape that exploits efficient and sophisticated techniques from the surface-integral-equation formulation 
of classical electromagnetic scattering. Unlike previous approaches to non-equilibrium fluctuations that involve 
scattering matrices — relating "incoming" and "outgoing" waves from each body — our approach is formulated 
in terms of "unknown" surface currents, laying at the surfaces of the bodies, that need not satisfy any wave 
equation. We show that our formulation can be applied as a spectral method to obtain fast-converging semi- 
analytical formulas in high-symmetry geometries using specialized spectral bases that conform to the surfaces 
of the bodies (e.g. Fourier series for planar bodies or spherical harmonics for spherical bodies), and can also 
be employed as a numerical method by exploiting the generality of surface meshes/grids to obtain results in 
more complicated geometries (e.g. interleaved bodies as well as bodies with sharp corners). In particular, our 
formalism allows direct application of the boundary-element method, a robust and powerful numerical imple- 
mentation of the surface-integral formulation of classical electromagnetism, which we use to obtain results in 
new geometries, such as the heat transfer between finite slabs, cylinders, and cones. 
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I. INTRODUCTION 

Quantum and thermal fluctuations of charges in other- 
wise neutral bodies lead to stochastic electromagnetic (EM) 
fields everywhere in space. In systems at equilibrium, these 
fluctuations give rise to Casimir forces (generalizations of 
van der Waals interactions between macroscopic bodies) that 
have recently become the subject of intense theoretical and 
experimental work.^~^ In non-equilibrium situations involving 
bodies at different temperatures, these fields also mediate en- 
ergy exchange from the hotter to the colder bodies, a process 
known as radiative heat transfer. Although the basic theoret- 
ical formalism for studying heat transfer was laid out decades 
ago"^"^ only recently have experiments reached the precision 
required to measure them at the microscale,^"^^ sparking re- 
newed interest in the study of these interactions in complex 
geometries that deviate from the simple parallel-plate struc- 
tures of the past.^^~^^ In this manuscript, we present a novel 
formulation of radiative heat transfer for arbitrary geometries 
based on the well-known surface-integral-equation (SIE) for- 
mulation of classical electromagnetism,^"^"^^ which extends 
our recently developed fluctuating surface-current (FSC) ap- 
proach to equilibrium Casimir forces^^ to the non-equilibrium 
problem of energy transfer between bodies of inequal tem- 
peratures. Unlike scattering formulations based on basis ex- 
pansions of the field unknowns best suited to special^^"^^ or 
non-interleaved periodic^^'^^"^^ geometries, or formulations 
based on expensive, brute-force time-domain simulations^^ 
and Green's functions calculations, "^^'"^^ this approach allows 
direct application of the boundary element method (BEM): 
a mature and sophisticated SIE formulation of the scattering 
problem in which the EM fields are determined by the solu- 
tion of an algebraic equation involving a smaller set of surface 
unknowns (fictitious surface currents in the surfaces of the ob- 
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A terse derivation of our FSC formulation for heat trans- 
fer was previously published in Ref. 42. The primary goals 
of this paper are to provide a more detailed presentation of 
this derivation and to generalize our previous formula for the 



heat transfer between two bodies to other situations of inter- 
est, including geometries consisting of multiple and/or nested 
bodies. We also demonstrate that the FSC framework can be 
applied as a spectral method to obtain semi-analytical for- 
mulas in special geometries with high symmetry, as well as 
for purely numerical evaluation using BEM, which we ex- 
ploit to obtain new results in a number of complicated ge- 
ometries that prove challenging for semi-analytical calcula- 
tions. Although our formulation here employs similar guiding 
principles as our previous work on equilibrium Casimir phe- 
nomena^^'^^ — both are based on the SIE framework of clas- 
sical EM scattering — the heat-transfer case is by no means a 
straightforward extension of force calculations, because gen- 
eralizing the equilibrium framework to non-equilibrium situ- 
ations requires very different theoretical techniques. For ex- 
ample, the fact that in Ref. 28 we considered only equilib- 
rium fluctuations made it possible for us to directly exploit 
the fluctuation-dissipation theorem for EM fields, "^"^ which re- 
lates the field-field correlation function at two points to a sin- 
gle Green's function between those two points. In contrast, 
although a fluctuation-dissipation theorem exists in the non- 
equilibrium problem, the field-field correlation functions are 
in this case determined by a product of two Green's functions 
integrated over the volumes of the bodies. ^^'"^"^ A key step in 
our derivation below is a correspondence between this vol- 
ume integral (involving products of fields) and an equivalent 
surface integral involving the fictitious surface currents and 
fields of the SIE framework, that was not required in the equi- 
librium case. 

The heat radiation and heat transfer of bodies with sizes 
and/or separations comparable to the thermal wavelength 
can deviate strongly from the predictions of the Stefan- 
Boltzmann law.^'^^ For instance, in the far field (object sep- 
arations d much greater than the thermal wavelength At = 
hc/kBT), radiative heat transfer is dominated by the ex- 
change of propagating waves and is thus nearly insensi- 
tive to changes in separations (oscillations from interfer- 
ence effects typically being small^"^^). In the less-studied 
near field regime (d < At), not only are interference ef- 



fects important, but otherwise-negligible evanescent waves 
also contribute flux through tunneling. ^^'^^ Such near-field 
effects have been most commonly studied in planar ge- 
ometries, where the monotonically increasing contribution 
of evanescent waves with decreasing d results in orders-of- 
magnitude enhancement of the net radiative heat transfer 
rate (exceeding the far-field black-body limit at sub-micron 
separations^^). This enhancement was predicted theoreti- 
cally^'^^'^^ and observed experimentally "^^""^^ decades ago in 
various planar structures, and has recently become the subject 
of increased attention due to its potential application in nan- 
otechnology, with ramifications for thermal photo voltaics^^'^^ 
and thermal rectification,^^"^^ nanolithography,^^ thermally 
assisted magnetic recording,^^ and high-resolution surface 
imaging. ^"^'^^'^^ Thus far, there have been numerous works fo- 
cused on the effects of material choice in planar bodies,^^'^^ 
including studies of graphene sheets,^^ hyperbolic^^ and 
anisotropic materials, ^^ and even materials exhibiting phase 
transitions,^^ to name a few. Along the same lines, many 
authors have explored transfer mediated by surface polari- 
tons in thin films^^"^^ and Id-periodic planar bodies. ^^ Despite 
decades of research, little is known about the near-field heat 
transfer characteristics of bodies whose shapes differ signif- 
icantly from these planar, unpatterned structures. Theoreti- 
cal calculations were only recently extended to handle more 
complicated geometries, including spheres,^^'^^ cylinders, ^^ 
and cones^^ suspended above slabs, dipoles interacting with 
other dipoles"^^'^^"^"^ or with surf aces, "^^'^^"^^ and also pat- 
terned/periodic surfaces . ^^'^ ^ 39,42,78-80 

General-purpose methods for modelling heat transfer be- 
tween bodies of arbitrary shapes can be distinguished in at 
least two ways, in the abstrsict formulation of the heat-transfer 
problem and in the basis used to "discretize" the formulation 
into a finite number of unknowns for solution on a computer 
(or by hand).^^ Theoretical work on heat transfer has mainly 
centered on "scattering-matrix" formulations which express 
the heat transfer in terms of the matrices relating incoming and 
outgoing wave solutions from each body.^^~^^'^^'^^'^^ These 
formulations tend to be closely associated with "spectral" dis- 
cretization techniques, in which a Fourier- like basis (Fourier 
series, spectral harmonics, etc.) is used to expand the un- 
knowns, because the incoming/outgoing waves must be ex- 
pressed in terms of known solutions of Maxwell's equations, 
which are typically a spectral basis of planewaves, spherical 
waves, and so on. Such a spectral basis has the advantage 
that it can be extremely efficient (exponentially convergent) if 
the basis is specially designed for the geometry at hand (e.g. 
spherical waves for spherical bodies^^). Scattering-matrix 
methods can also be used for arbitrary geometries, e.g. by ex- 
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or 



panding arbitrary periodic structures in Fourier series 
by coupling to a generic grid/mesh discretization to solve the 
scattering problems, ^^'"^^'^^ but exponential convergence no 
longer generally obtains. Furthermore, Fourier or spherical- 
harmonic bases of incoming/outgoing waves correspond to 
uniform angular/spatial resolution and require a separating 
plane/sphere between bodies, which can be a disadvantage 
for interleaved bodies or bodies with comers or other fea- 
tures favoring nonuniform resolution. In contrast to the geo- 



metric specificity encoded in a particular scattering basis, one 
extremely generic approach is a brute-force discretization of 
space and time, allowing one to solve for heat-transfer by a 
Langevin approach^^ that handles all geometries equally, in- 
cluding geometries with continuously varying material prop- 
erties. The FSC approach lies midway between these two 
extremes. Like the scattering-matrix approach, the FSC ap- 
proach exploits the fact that one knows the EM solutions 
(Green's functions) analytically in homogeneous regions, so 
for piecewise-homogeneous geometries the only remaining 
task is to match boundary conditions at interfaces. Unlike 
the scattering-matrix approach, however, the FSC approach is 
formulated in terms of unknown surface currents rather than 
incoming/outgoing waves — the surface currents are arbitrary 
vector fields and need not satisfy any wave equation, which 
leads to great flexibility in the choice of basis. As described 
in this paper, the FSC formulation can use either a spectral 
basis or a generic grid/mesh and, as demonstrated in Ref. 42 
and Ref. 80, works equally well for interleaved bodies (lack- 
ing a separating plane or even a well-defined notion of "in- 
coming/outgoing" wave solutions). Moreover, the FSC for- 
mulation reduces the heat-transfer problem to a simple trace 
formula in terms of well- studied matrices that arise in SIE for- 
mulations of classical EM, which allows mature BEM solvers 
to be exploited with minimal additional computational effort. 
The radiative heat transfer between two bodies 1 and 2 at 
local temperatures T^ and T^ can be written as:^^'^^ 



H= duj [e(cj, T^) - e(cj, T^)] <^( 
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(1) 



where 0(co',T) = huj/[ex.-p{huj/kBT) — 1] is the Planck en- 
ergy per oscillator at temperature T, and ^ is an ensemble- 
averaged flux spectrum into body 2 due to random currents 
in body 1 (defined more precisely below via the fluctuation- 
dissipation theorem"^'"^^'^^). (Physically, there are currents in 
both bodies, but EM reciprocity^^ means that one obtains the 
same ^ for flux into body 1 from sources in body 2; this 
also ensures that H obeys the second law of thermodynam- 
ics.) The only question is how to compute ^, which naively 
involves a cumbersome number of scattering calculations. 

The main result of this manuscript is the compact trace- 
formula for ^ derived in Sec. II, which involves standard ma- 
trices that arise in BEM calculations and forgoes any need 
for evaluation of fields or sources in the volumes of the bod- 
ies, separation of incoming and outgoing waves, integration 
of Poynting fluxes, or many scattering calculations. As ex- 
plained below in Sec. HID and Sec. IIIC, by a slight modifi- 
cation of the two-body formula one can also straightforwardly 
compute the spatially resolved pattern of Poynting flux on the 
surfaces of the bodies, as well as the emissivity of an isolated 
body. Section III A illustrates how important physical proper- 
ties such as reciprocity and positivity of heat transfer manifest 
in the algebraic structure of the formulas. In Sec. HIE, we 
generalize the two-body formula to also describe situations 
involving multiple and/or nested bodies. The remaining sec- 
tions of the paper are devoted to validating the FSC formal- 
ism by checking it against known results in special geome- 
tries consisting of spheres and semi-infinite plates, as well as 



applying it to obtain new results in more complicated geome- 
tries consisting of finite slabs, cylinders, and cones. Specif- 
ically, Sec. IV B considers application of the FSC formula- 
tion in high-symmetry geometries where the use of special- 
bases expansions involving Fourier and spherical-wave eigen- 
functions (provided in Appendix A) leads to fast-converging 
semi-analytical formulas of heat radiation and heat transfer 
for spheres and semi-infinite plates. In Sec. IV C and Sec. V, 
we exploit a sophisticated numerical implementation of the 
FSC formulation based on BEM to check the predictions of 
the semi-analytical formulas in the case of spheres and to ob- 
tain new results in more complex geometries. Finally, the ap- 
pendices at the end of the paper provide additional discussions 
that supplement and aid our derivations in Sec. II and Sec. III. 
Specifically, Appendix B provides a concise derivation of the 
principle of equivalence and its application to SIEs, and Ap- 
pendices C 1 and C 2 provide succinct proofs of reciprocity 
and positivity of Green's functions and SIE matrices, respec- 
tively. 



II. FSC FORMULATION 

In this section, we review the SIE method of EM scatter- 
ing and apply it to derive an FSC formulation of radiative heat 
transfer between two bodies. The result of this derivation is 
a compact trace expression for ^ involving SIE matrices. We 
further elaborate on these results in Sec. Ill, where we extend 
the formulation to handle other situations of interest, includ- 
ing the emissivity of isolated bodies, distribution of Poynting 
flux on the surfaces of the bodies, and heat transfer between 
multiple and/or nested bodies. 



A. Notation 



Let 6 = 



denote 6-component 



h) ^"'i^ = (k 

volume electric and magnetic fields and currents, respectively, 
and ^ denote 6-component surface currents (which technically 
have only 4 degrees of freedom since they are constrained to 
flow tangentially to the surfaces). In a homogeneous medium, 
fields are related to currents via convolutions (^) with a 6 x 6 
homogeneous Green's tensor r(x, y) = r(x — y,0), such 
that = r^((T + ^),or more explicitly 
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dVr(x,y)[a(y)+^(y)], 
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where 






ik 



-C i^ 



is the Green's tensor composed of 3 x 3 electric and magnetic 
Dyadic Green's functions (DGFs), determined by the "pho- 
ton" DGFs G and C. In the specific case of isotropic media 
(scalar e and /i), G and C satisfy 




Figure 1. Schematic depicting two disconnected bodies described by 
surfaces dV^ and dV^ , and held at temperature T^ and T^, respec- 
tively. Surface currents ^^ and ^^ laying on the surfaces of the bodies 
give rise to scattered fields (j)^~ and 0^", respectively, in the interior 
of the bodies, and scattered field (jp~ in the intervening medium 0. 



and C = |V X G, with wavenumber k = ujy/e]l and 
impedance Z = \J ^je. Our derivation below applies to ar- 
bitrary linear anisotropic permittivity e and permeability /i, 
so long as they are complex- symmetric matrices in order to 
satisfy reciprocity^"^ (see Appendix CI). The mathematical 
consequence of reciprocity, as described in the Appendix, is 
that F is complex-symmetric up to sign flips. In particular, 
r(x,x')^ = 5r(x,x')5, where the 6 x 6 matrix S = 5"^ 
flips the sign of the magnetic components. This reciprocity 
property is a key element of our derivation below. 



B. Surface integral equations 

Consider the system depicted in Fig. 1, consisting of two 
homogeneous bodies, 1 and 2 (volumes V^ and V'^ and tem- 
peratures T^ and T^), separated by a lossless medium (vol- 
ume V^) by two interfaces dV^ and dV^, respectively. Con- 
sider also sources a^ located in the interior of V^, and de- 
note the total fields in each region by (j)^. The homogeneous- 
medium Green's functions for the infinite media in region r 
are denoted by F^. Consider also the decomposition of the 
total fields 0^ in each region r into "incident" fields ^^+ (due 
to sources within r) and "scattered" fields (j)^~ (from interac- 
tions with the other regions, including both scattering off the 
interface and sources in the other regions). That is, we can 
write: 0^ = ^^+ + ^^-, with 0^+ = F^ ^ a^ 

The core idea in the SIE formulation is the principle of 
equivalence, ^^'^^~^^ whose derivation is briefly reprized in Ap- 
pendix B, which states that the scattered field ^^~ can be ex- 
pressed as the field of some fictitious electric and magnetic 
surface currents ^^ located on the boundary of region r, act- 
ing within an infinite homogeneous medium r. In particular, 
one can write: 






(4) 
(5) 



[V X V X - fc2] G(fc; X, x') = S{x - x')I, 



(3) 



for r = 1,2, with fictitious currents ^^ completely determined 
by the boundary condition of continuous tangential fields at 
the body interfaces. Specifically, equating the tangential com- 
ponents of the total fields at the surfaces of the bodies, we find 



the integral equations: 



A^+ 
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lay^ 



IdV^ 



(6) 



which can be solved to obtain ^^ from the incident fields. 
This is the "PMCHW" surface-integral formulation of EM 
scattering. ^"^'^^'^^ 

Let {/3^} be a basis of 6-component tangential vector fields 
on the surface of body r, so that any surface current ^^ can be 
written in the form ^^(x) = ^^ x5^/3^(x) for N coefficients 
{x^}. (In BEM, f3n is typically a piecewise-polynomial "ele- 
ment" function defined within discretized patches of each sur- 
face, most commonly the "RWG" basis functions^^.^^ How- 
ever, one could just as easily choose f3n to be a spherical har- 
monic or some other "spectral" Fourier-like basis, as shown 
in Sec. IV B. The key point is that (3^ is an arbitrary basis of 
surface vector fields; unlike scattering-matrix formulations, it 
need not consist of "incoming" or "outgoing" waves nor sat- 
isfy any wave equation.) Taking the inner product of both 
sides of Eq. 6 with (3^ (a Galerkin discretization^"^), one ob- 
tains a matrix "BEM" equation of the form: 
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where x 



represents the expansion of the surface 



currents ^^ = En</^n"5 



the incident fields sL 
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describes interactions with matrix elements G^^ = (/3^, F^^ 
/3l) among basis functions. G^ represents multi-body interac- 
tions between basis functions on both bodies, via waves prop- 
agating through the intervening medium 0. G^ represent self- 
interactions via waves propagating within a body, given by: 



Gr 



/^r,rr /or j^r . R'^\ 



(9) 



Here, (•,•) denotes the standard inner product {(p^tl^) = 
f (p'^^l^, with the * superscripts denoting the conjugate- 
transpose (adjoint) operation. 

A key property of the Green's function is reciprocity, as 
summarized and derived in Appendix C 2, and this property is 
reflected in symmetries of the matrices G and W. For sim- 
plicity, let us begin by considering the case of real- valued 
basis functions /3n- Let S be the matrix such that Sx flips 
the signs of the magnetic components (assuming that we ei- 
ther have separate basis functions for electric and magnetic 
components, as in the RWG basis, or more generally that 
the basis functions come in f3n and 8/3^ pairs). Note that 
S~^ = S = S^ = S*. In this case, as reviewed in Ap- 
pendix C2, it follows that W^ = SWS and G^ = SGS. 



Once we have derived our heat-transfer formula for such real- 
valued basis functions, it is straightforward to generalize to 
complex- valued bases as described in Sec. IIIB. 



C. Flux spectrum 

Our goal is to compute the flux spectrum <l> into V'^ (the 
absorbed power in body 2) due to dipole current sources a^ 
in V^ (integrated over all possible positions and orientations). 
We begin by considering ^^n, or the flux into body 2 due 
to a single dipole source a^ within body 1, corresponding to 
(^^+ = F^ ^ a\ with 00+ = 02+ = 0. In the SIE Eq. 7, 
this results in a source term s with s^ = (/3^,F^ ^ a^) 
and s^ = 0. As derived in Appendix B, the Poynting flux 
can be computed using the fact that ^ is actually equal to the 

surface-tangential fields,^^ ^ 



^ , , where n is the 
^-n X Ey 

outward unit-normal vector. It follows that the integrated flux 
-^ Re ^(E X H) • n = I Re(^^ 0°). (This can also be de- 
rived as the power exerted on the surface currents by the total 
field, with an additional 1/2 factor from a subtlety of evaluat- 
ing the fields exactly on the surface. ^^) Hence: 



^n 



1 



\Re{e.^') = -^Re{e 



where we used the continuity of 
02+ = 0. Substituting ^^ = ^^ , 
nition of G^ in Eq. 8, we obtain 



02) = iRe{e^-^2*a, 

0^ and 02 and the fact that 
v'^Pn ^^^ recalling the defi- 



^n 



1 



Re (x2*G2x2) = - J Re (x'^G^x) 

1 
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1 
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(symG'^) X = -jS*W* UymG'^) Ws 



Tr 



ss'^W 



(symG^) 



W 



where sym G = ^{G -{-G"^) denotes the Hermitian part of G. 
Computing ^^i is therefore straightforward for a single 
source a^. However, the total spectrum 



^ = {^^i) = -^ Tr {ss*)W* (symG^^ W 



(10) 



involves an ensemble-average (• • • ) over all sources a^ and 
polarizations in F^. While this integration can be performed 
explicitly, we instead seek to simplify matters so that the fi- 
nal expression for $ involves only surface integrals. The key 
point is that ss"^ is sm N x N matrix describing interactions 
among the N surface-current basis functions. The ensemble 
average (ss*) is also an A^ x TV matrix, which we would like to 
express in terms of a simple scattering problem involving the 
SIE Green's function matrices, hence eliminating any explicit 
computations over the interior volume V'^. 

Defining the Hermitian matrix C = (s5*), it follows 
that its only non-zero entries lie in the upper-left Ni x A^i 
block C^ = (5^5^*), and are given by C^^ = {sl^s}^) = 
{{Pi^T^ ^ a') {T^^ct\ Pi)), or 
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d^^JJJ cPyPl{^rT\^,y)a\y)(jj)cP^ jjj dVV(y')*ri(x',y')*/3i(x') 



(11) 



where in the third Hne we have performed an integration over 
all dipole positions by employing the fluctuation-dissipation 
theorem^^ for the current-current correlation function, 

(^'(y)^'(y')*) = -L0lmxiy,uj)5iy-y'), (12) 

TT 

and where we omitted the dependence on the Planck energy 

distribution 6(co', T), which has been factored out into Eq. 1, 

and where Im x denotes the imaginary part of the 6 x 6 ma- 

/im e 
terial susceptibility tensor, so that Imx = ( q t 

which is related to material absorption. 



Equation 11 closely resembles an absorbed power in the 
volume of body 1, since absorbed power for a field ^ is 
^ / 0* (cj Im x)0-^^ To make this analogy precise, some care- 
ful algebraic manipulation is required, and the abovemen- 
tioned reciprocity relations [r(x, x')^ = «Sr(x, x')5, W^ = 
SWS, etc] play a key role. In particular, the fact that C^ is 
Hermitian implies that the matrix is completely determined 
by the values of x^*5'(C^)^5'x^ for all x^, where we have in- 
serted the sign-flip matrices S and the transposition for later 
convenience. Interpreting x^ as the basis coefficients of a sur- 
face current ^^ = ^^ x^^Pn ^^ dV^, we find: 



x^^'SiC^fSx^ = 



x^'^Ss^ 



(jj)d'^jjjd'yjj)d'^jjjd'y'e{^rSf^^ 

l§d^^jjj d'y jj) d^^e (x) *5fn^^ [^ Im X(y)] SV' (x^ y)e (x^ 



- # d^x jjj d'y (jj) d^^e (x) *5ri (x, y)S [uj Im x(y )] T' (x^ y)e (x^ 
d^x fffd'ySd'^' [ri(y,x)ei(x)]* [c^Imx(y)] [T\^\y)e{^')] 



-(rUe\(^imx)r^^^^) 



(13) 



where in the first and fourth lines we invoked reciprocity 
(from above) and in the third line we assumed that S com- 
mutes with Imx, which is true for reciprocal media. (The 
only way that S would not commute with Im x would be if 
there were a chiral susceptibility coupling electric and mag- 
netic fields directly, also called a bi-anisotropic susceptibility, 
which breaks reciprocity. ^^'^^) Letting (f)^ =r^^^^be the 
field due to the surface current ^^, it follows that 



x^^'SiC^fSx^ 
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b\{ojlmx)^'). 



(14) 



But, as noted above, ^ 



^ , cj (Im x) 0^ ) (where the inner prod- 
uct (•,•) is now over the volume V^) has a simple meaning: it 
is the absorbed power in V'^ from the currents ^^, or equiva- 
lently, the time-average power density dissipated in the inte- 



rior of body 1 by the field (j)^ produced by ^^ . 

Computing the interior dissipated power from an arbi- 
trary surface current turns out to be somewhat complicated, 
since one needs to take into account the possibility that 
the equivalent surface currents arise from sources both out- 
side and inside V'^. If, on the other hand, we could re- 
strict ourselves to equivalent currents ^^ which are out- 
side of F^, then we can use the result from above that 
the incoming Poynting flux (the absorbed power) is sim- 



ply-lRe{e.^') 



i\ - 



1^1* 



(symG^)x^. Substituting this 



into Eq. 14, we would be immediately led to the identity 
x^''S{C^)^Sx^ = -f Re (xi*G^x^), and this gives an ex- 
pression for C^ in terms of G^. It turns out that indeed, we 
need not handle arbitrary ^^ since the C matrix is never used 



by itself — it is only used in the trace expression 



<!> 



Tr[ 



SWSS{sym G^)SSW*SC'^ 
Tr[SW{symG^)W*SC^] 

Tr[SC'^ SW{symG^)W*], 



(15) 



using reciprocity. As shown in Sec. Ill A, the standard defi- 
niteness properties of the Green's functions (currents do non- 
negative work) imply that sym G^ is negative semidefinite and 
hence admits a Cholesky factorization^^ symG^ = —tj^*tj^. 
It follows that Eq. 15 can be written as -\ Tr[X*S'C^S'X], 
where X = Wtj^* are the "currents" due to "sources" rep- 
resented by the columns of t/^*, which are all of the form 

: currents from sources in V'^ alone. So, effectively, 



SC^S is only used to evaluate the power dissipated in V^ 
from sources in V'^, and by the same Poynting-theorem rea- 
soning from above, it follows that S{C^)^S = — f symG^, 
and hence 



C 



■ sym S{G^)'S 



sym G^ 



(16) 



by the symmetry of &. Substituting this result into Eq. 10 
then gives the heat transfer formulation summarized in the 
next section. 



D. Heat transfer formula 

The result of the above derivation is that the ensemble- 
averaged flux from F^ to F^ can be expressed in the compact 
form: 



27r 



(symG^W* (symG^W 



27r 



Tr [(symG^) W^^* (symG^) W^^] 



(17) 
(18) 



with W^^ relating incident fields at the surface of body 2 to 
the equivalent currents at the surface of body 1 . Our simpli- 
fied expression is computationally convenient because it only 
involves standard matrices that arise in BEM calculations,^^ 
with no explicit need for evaluation of fields or sources 
in the volumes, ^^'^^"^^ separation of incoming and outgoing 
waves, ^^~^^'^^'^^'^^ integration of Poynting fluxes, ^^ or any ad- 
ditional scattering calculations. 



III. GENERALIZATIONS 

In this section, we study the positivity and symmetries of 
the two-body heat-transfer formula above and consider gen- 
eralizations to include other situations of interest. Following 



similar arguments as those employed in the previous section, 
we derive formulas for the emissivity of isolated bodies, the 
spatial distribution of Poynting flux on the surfaces of bod- 
ies, and the heat transfer between multiple and nested bodies. 
In Sec. IIIB, we show that abandoning our choice of real-/3 
basis functions above in favor of complex-/^ functions does 
not change the final formula for ^, so long as the /3s come in 
complex conjugate pairs. 



A. Positivity and reciprocity 

In addition to its computational elegance, Eq. 18 alge- 
braically captures crucial physical properties of the flux spec- 
trum: ^ is positive-definite ^ > and symmetric with respect 
to 1 o 2 exchange, as required by reciprocity. Of course, the 
positivity of ^ is immediately clear from the Rytov starting 
point of fluctuating currents inside the bodies: the absorbed 
power in one body from sources in the other body is simply 
~ /(a;Ime)|Ep > (since ujlme > for passive me- 
^j^83,84^ Hence, positivity must hold for any formulation that 
is mathematically equivalent to the Rytov picture. However, 
it is still useful and nontrivial to understand how this positiv- 
ity manifests itself algebraically in a given formulation. For 
example, Ref. 35 showed how positivity manifests itself in a 
scattering-matrix framework. In our ESC framework, positiv- 
ity turns out to correspond to the fact that $ can be interpreted 
as kind of matrix norm. 

As derived above, the standard definiteness properties of 
the Green's functions (currents do nonnegative work) im- 
ply that symG^ is negative semidefinite and hence admits 
a Cholesky factorization symG^ = —U^*U^, where U^ is 
upper-triangular. It follows that 

<!> = ^ Tr \U^WU'^*U'^WU^''] 
27r ^ ^ 



^Tr[Z*Z] = ^||Z||^, 



(19) 



where Z = U^WU^'^, is a weighted Erobenius norm of the 
SIE matrix W, which from above we know is necessarily non- 
negative. 

Furthermore, reciprocity (symmetry of ^ under 1 o 2 in- 
terchange) corresponds to simple symmetries of the matri- 
ces. As derived in Appendix CI, r(y,x)^ = 5r(x, y)5, 
G^ = SGS and W^ = SWS, where S = S^ = S'^ = S"" 
is the matrix that flips the signs of the magnetic basis coef- 
ficients and swaps the coefficients of /3n and /3n- It follows 
that 



<!> = ^Tr 

27r 



27r 



Tr 



SWS (sym SG^S^ SWS (sym SG^S^ 
TsymG^W* (symG^)w 



(20) 



where the 5* factors cancel, leading to the 1 ^ 2 exchange. 



B. Complex- valued basis functions 

For convenience, we assumed above that the basis functions 
Pn were purely real- valued. However, it easy to generalize 
the final result a posteriori to complex- valued basis functions. 
The relevant case to consider are basis functions that come in 
complex-conjugate pairs /3n and Pn' = Pn (true for any prac- 
tical complex basis). Such a basis can always be transformed 
into an equivalent real-valued basis Pn by the linear transfor- 
mation pn = 7i(/^n + Pn') and pn' = ^if^n - Pn')- In an 

expansion ^ = Y.n ^nPn = Zln ^nPn. this is simply a rota- 
tion X = Qx where the matrix Q is easily verified to be uni- 
tary (Q* = Q~^), since it is composed of unitary 2x2 blocks 
(operating on n^n' complex-conjugate pairs). Given such a 
unitary change of basis, we can make a corresponding unitary 
change to the G and W matrices from above. G = QGQ* 
and W = QWQ*, to obtain the matrices in the complex ba- 
sis. By inspection of the ^ expression above, all of the Q 
factors cancel after the change of basis and one obtains the 
same expression in the complex basis with the new G and W 
matrices. 



C. Emissivity of a single body 

The same formalism can be applied to compute the emis- 
sivity of a single body. For a single body 1 in medium 0, the 
emissivity of the body is the flux <l>^ of random sources in V'^ 
into V^ P Following the derivation above, the flux into V^ is 



4Re(e\( 



4 (^\ r° ^ ^^ ) . The rest of the derivation 



is essentially unchanged except that W = {G^ -\-G^)~ 
there is no second surface. Hence, we obtain 



1 



<^° = -- Tr [(symG^) W* (symG°) W] , 



27r 



smce 



(21) 




Figure 2. Schematic depicting three disconnected bodies described 
by surfaces dV^, dV'^, and dV^, and held at temperature T\ T^ 
and T^, respectively. Surface currents ^^, ^^, and ^^, laying on the 
surfaces of the bodies give rise to scattered fields 0^" , 0^" , and 0^" , 
respectively, in the interior of the bodies, and scattered field 0°" in 
the intervening medium 0. 



Note that $ = TrF^. Similarly, by swapping 1 ^ 2 we 
obtain a matrix F^ such that ^n = F^^ is the contribution of 
/3^ to the flux on surface dV^. 



E. Multiple and nested bodies 

In this section, we extend the FSC formalism above to sit- 
uations involving multiple and nested bodies. For simplicity, 
we only consider an additional medium 3, since generaliza- 
tions to include additional bodies or levels of nesting readily 
follow. Because the derivation is almost identical to the two- 
body case, we only focus on those aspects that differ. 



which again is invariant under 1^0 interchange from the 
reciprocity relations (Kirchhoff 's law). 



D. Surface Poynting-flux pattern 

It is also interesting to consider the spatial distribution of 
Poynting-flux pattern, which can be obtained easily because, 
as explained above, | Re[^^(x)*^^(x)] is exactly the inward 
Poynting flux at a point x on surface 2. It follows that the 
mean contribution <l>^ of a basis function /3^ to ^ is 



^^ = 



4 

= ^Re 

27r 



Re s'^W'^elel'^G^Ws 
el''G^wUymG^)Wel 



where e^ is the unit vector corresponding to the /3^ compo- 
nent. This further simplifies to <l>^ = F^^, where 



^Re[G^W^^symG^)W^^''] . 



(22) 



1. Multiple bodies 

Consider the system depicted in Fig. 2, consisting of three 
disconnected bodies at different temperatures. Applying the 
principle of equivalence, one finds: 

(/>« = (/>«++ F^ ^(^1+ ^2 +f), 

for r = 1,2,3, with fictitious currents ^^ determined by the 
boundary conditions of continuous tangential fields at the the 
body interfaces. Equating the tangential components of the 
fields at the surfaces of the bodies, one obtains the integral 
equations: 



(rVr^)^e^+y(ro^e) 






^r+ 



aO+I 



IdV^ 



IdV^ 



(23) 




Figure 3. Similar three-body geometry as that depicted in Fig. 2, 
but with body 3 now embedded in the interior of body 2. Here, the 
scattered field (j)^~ includes contributions from both ^^ and ^^. 



along with the corresponding SIE matrix: 



V^ll VI^12 |yl3 
^21 1^22 1^23 
lySl |y32 1^33 



w-i 



^0,11 (^0,12 (^0,13 

_ I (70,21 (^0,22 (^0,23 

(^0,31 (^0,32 (^0,33 
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G^ 



G^ 



I . (24) 

G^ 



Gl (72 (73 

The derivation of the flux spectrum for any given pair of bod- 
ies mirrors exactly the derivation in Sec. II, with the only dif- 
ference being the modified SIE matrix W . The final expres- 
sion for the flux spectrum into V^ due to random currents in 
y*/j is given by: 

^'^ = — Tr [(symff )W^'^*(symG^')iy^'^] , (25) 

which again is invariant under i ^ j interchange. 



2, Nested bodies 

Consider now the system depicted in Fig. 3, involving 
three bodies at different temperatures with one of the bod- 
ies (medium 2) containing another (medium 3). Applying the 
principle of equivalence again, one finds: 



^^ = ^^+-r^^(^^-a, 

for r = 1,3, with fictitious currents ^^ determined by the 
boundary conditions of continuous tangential fields at the 
body interfaces. Equating the tangential components of the 
fields at the the surfaces of the bodies, one obtains the integral 
equations: 

- ^2+ _ / 0+ I 






\dV^ 



I ay 2 5 



(r^ + r^)^^ -r^^e^ 



a2+| 



where dV^ denotes the interface between V'^ and V^ , from 
which one obtains the corresponding SIE matrix: 



lyll VI^12 y^l3 
1^21 1^22 ^23 
1^31 1^32 1^33 



(^0,11 (70,12 

(70,21 (70,22 



w- 



GO 



G^ 



G^ 

_(72,32 



_^2,23 
(72,33 



G^ 



Gi 



G2 



G3 



(26) 



Although the derivation of the flux spectrum for any given 
pair of bodies closely mirrors the derivation in Sec. II C, im- 
portant deviations arise due to the difference in topology. In 
what follows, we only focus on those steps that differ signif- 
icantly. The asymmetry of the geometry also means that we 
must consider $ for each pair separately. 

First, we compute the flux spectrum $^^ into V^ (the ab- 
sorbed power in 3) due to dipole current sources in F^. The 
flux into body 3 due to a single dipole source cr^ inside body 
1 is given by: 



$ 



13 



1 



Ke(r,r) = 7l^e(^ ' 



Re{e 



1 



_r3- 



3* /o3^3\ 



Re [x^'^G'^x 



After ensemble averaging over a^ as before, we obtain: 
1 



^ 



13 



27r 



Tr [(symG^)W^^*(symG^)iy' 



/311 



(27) 



Second, we compute the flux spectrum <l>^^ into V'^ (the ab- 
sorbed power in body 2) due to dipole current sources in F^. 
Direct application of Poynting's theorem at dV'^ in this case 
does not yield <l>^^ but rather the quantity we denote as <l>^(^) : 
the flux into the entire region contained by dV^ from sources 
in V^, which includes absorption in both V'^ and V^ . It fol- 
lows that ^^^ = <l>^^^) — <l>^3. So, it only remains to compute 
(^^^'^\ starting with the flux from a single cr^ source, given by: 



= ^Ree 



-r^ * {^ 



-iRe[.-(G^ 



2^2 






with the additional x^ term stemming from absorbed power in 
body 3. We ensemble average as before, and obtain 



^ 



12 



<|)1(2) 
1 



$ 



13 



27r 



Tr 



(symG^)iy^^*(symG^)iy 



21 



\dV^ 



\dV^ 



(sym G^) sym (iy2^*G^'^^ W' 



31^ 



^ 



13 



(28) 



Finally, we compute the flux spectrum <^^^ into V^ (the ab- 
sorbed power in body 2) due to dipole current sources in V^. 
^^^ can be computed by subtracting the flux <l>^'^^) leaving 
body 2 through dV'^ from the flux ^^'^^^ entering body 2 
through dV^. Specifically, for a single dipole a^, we find: 

= \ Re{e, -T' H')-\ Re(e^ -T^ * {f - f)) 



Re {x^*G^x^) + ^ Re [x^* {G^x^ - G^'^^x"")] . 



The final result is the expression: 

^32 ^^3(3) _^3(2)^ ^29) 

with 

$3(2) ^ ^xr r(symG3)V[/23*(symG2)W23 

27r 

+ (symG^)sym(W^^*G^'2^W^^)] (30) 

^^(^^ = ^Tr [(symG^)W^^*(symG^)W^^] . (31) 

27r 

For example, the heat transfer between V'^ and the com- 
bined V^'^'> = F^ U F^ is given by 

H^^^^ = /Gti^'^'^ - Bt^^^' - eT3^^^ (32) 

where the integral is taken over all positive frequencies uj and 
6t = &{uJjT). In the special case T^ = T^, the expression 
reduces to the expected form: 



H'^^^= Aoti - OtO ^'^'^ 



(33) 



As before, we obtain reciprocity relations $*^ = ^^* between 
every pair of bodies, but these relations are no longer apparent 
merely by inspection of ^*^ . Because each body is topolog- 
ically distinct, ^^^ is no longer obtained from ^*^ merely by 
interchanging i and j, but instead must be derived separately 
(using analogous steps). Upon carrying out this derivation, 
we verify that $*-^ = ^^^ as required. Furthermore, the pos- 
itivity of $*-^ appears harder to derive algebraically from the 
final expression than in the non-nested cases, and we do not 
do so in this work. (Although it follows from the second law 
of thermodynamics, the scattering-matrix proof of positivity^^ 
should apply to nested bodies with minimal modification). 



IV. VALIDATION 

We now apply our FSC formulation to obtain results ob- 
tained previously using other scattering formulations in sev- 
eral high-symmetry geometries. In Sec. IV A, we discuss 



the choice of basis, contrasting BEMs that use a generic sur- 
face mesh with spectral methods that use a Fourier-like ba- 
sis, and point out that the latter are actually closely related to 
scattering-matrix methods in the case of high- symmetry ge- 
ometries. In Sec. IVB, we derive semi-analytical expressions 
of heat radiation and heat transfer for spheres and plates, us- 
ing surface spherical-harmonics and Fourier bases to describe 
the SIE surface unknowns, and show that these agree with 
previous formulae derived using other formulations. ^'^^'^^'^^ 
In Sec. IV C, we present a general-purpose numerical im- 
plementation of the FSC formulation based on a standard 
triangular-mesh discretization of the surfaces of the bodies 
known as the BEM "RWG" method; we check it against pre- 
vious heat-transfer methods by computing the heat transfer 
between spheres. 



A. Choice of basis 

The standard approach for solving the SIEs above is to dis- 
cretize them by introducing a finite set of basis functions /^^ 
defined on the surfaces of the bodies. As noted above, an 
important property of SIE formulations is that Pn is an arbi- 
trary basis of surface vector fields: unlike scattering-matrix 
formulations, ^^~^^ they need not satisfy any wave equation, 
nor encapsulate any global information about the scattering 
geometry, nor consist of "incoming" or "outgoing" waves into 
or out of the bodies. This lack of restriction on j3n is a power- 
ful property of the SIE formalism. 

There are two main categories of basis functions that one 
could employ: spectral bases or boundary -element bases. A 
spectral basis consists of a Fourier-like complete basis of non- 
localized functions, such as spherical harmonics or Cheby- 
shev polynomials,^"^ which are truncated to obtain a finite ba- 
sis. BEMs instead first discretize each surface into a mesh 
of polygonal elements (e.g. triangles), and describe functions 
piecewise by low-degree polynomials in each element.^"^'^^'^^ 
Spectral bases have the advantage that they can converge ex- 
ponentially fast for smooth functions, ^^ or in this case for 
smooth interfaces, but they are not as well suited to han- 
dle singularities such as comers, and moreover represent sur- 
faces with essentially uniform spatial resolution. A BEM ba- 
sis, on the other hand, is more flexible because it can use a 
nonuniform mesh to concentrate spatial resolution where it is 
needed,^^'^^ and furthermore the localized nature of the basis 
functions has numerical advantages in assembling and apply- 
ing the W and G (Green's function) matrices. ^^'^^^ The most 
common BEM technique employs a mesh of triangular ele- 
ments (panels) with vector- valued polynomial basis functions 
called an RWG (Rao-Wilson-Glisson) basis, ^^ where each 
basis function is associated with each edge of the mesh and 
is nonzero over a pair of triangles sharing that edge. Many 
years of research have been devoted to the efficient assembly 
of the G matrices for the RWG basis (by evaluating the singu- 
lar panel integrals of p),^^^"^^^ and to fast methods for solving 
the resulting linear equations. ^^^'^^"^ 

For a handful of highly symmetric geometries, however, 
spectral bases have an additional advantage: a special basis 
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can be chosen such that most of the matrix elements can be 
computed analytically (and many of the G matrices are di- 
agonal as a consequence of orthogonality). This has a close 
connection to scattering-methods, because whenever there is a 
known incoming/outgoing wave basis (e.g. spherical waves), 
one can construct an equivalent set of surface-current basis 
functions (e.g. spherical harmonics) by the principle of equiv- 
alence. (In fact, the principle of equivalence can be used to 
derive an exact equivalence between our ^ expressions and 
the analogous expressions from the scattering-matrix formu- 
lation, which we do not show here.) In the example of inter- 
actions between two spherical bodies, if we employ a (vec- 
tor) spherical-harmonic basis on each body, then the G^ self- 
interaction matrices are diagonal and the G^'^^ interaction 
matrix is given by "translation matrices" that relate spherical- 
wave bases at different origins. ^^^ In this way, by choosing a 
geometry- specific basis, the FSC formulation can retain all of 
the efficiency of the scattering-matrix methods, while preserv- 
ing the flexibility to employ a different basis as needed. 



B. Spectral basis 

In this section, we explicitly apply our FSC formulation 
with a spectral basis in three high- symmetry geometries for 
which the matrix elements can be evaluated semi-analytically: 
radiation of an isolated plate and an isolated sphere, and heat 
transfer between two parallel plates. In each case, we re- 
produce known solutions that were derived previously using 
scattering-matrix formulations. ^'^"^'^^'^^ The main purpose of 
this section is to illustrate how the FSC formulation with a 
spectral basis allows semi-analytical calculations similar to 
scattering-matrix formulations (albeit only in the handful of 
high- symmetry geometries where exact wave solutions can be 
constructed in each body). To begin with, we review the well- 
known spectral representation of the homogeneous DGF F in 
bases specialized to particular coordinate systems. 



1. Basis of Helmholtz solutions 

We wish to work with solutions of Maxwell's equations 
known analytically within each body and which are orthog- 

I 



onal when evaluated on the interfaces. These solutions, eval- 
uated at the interface of each body, will then provide a ba- 
sis of surface-tangential vector fields in which the G matrices 
can be evaluated analytically or semi-analytically. In partic- 
ular, we wish to work with solutions M and N of the vector 
Helmholtz equation (equivalent to Maxwell's equations in a 
homogeneous isotropic medium), ^^^ 



M 

N 



0, 



(34) 



with M = -i/k\/ X N and N = i/k\/ x M denoting 
purely electric and purely magnetic vector fields. (Note that 
M and N come in two flavors, depending on whether on 
solves Eq. 34 for outgoing or incoming boundary conditions.) 
Furthermore, we seek solutions of Eq. 34 in a coordinate sys- 
tem that allows separation of variables into "normal" and "tan- 
gential" components to some surface dV (which is possible 
for a small number of coordinate systems). We let r]± rep- 
resent the separable coordinate identified as the normal co- 
ordinate, and let r^ii represent the remaining tangential coor- 
dinates. The choice of coordinate system ultimately corre- 
sponds to a choice of basis, or independent solutions labeled 
by an index n that correspond to different scattering channels. 
Specifically, one is led to vector fields :^^^ 

Mt{v±.V\\) = <^(r?±)Xn(r7||) (35) 

with K,^ and {X^,Z^} denoting the normal and tangential 
components of the fields, and with ± denoting incoming (+) 
and outgoing (— ) solutions. For example, solutions in spher- 
ical coordinates yield the well-known vector spherical-wave 
solutions M^^(r, 6>,0) = i?^(r)Y^^^(6>, 0), described by 

spherical Hankel functions tzf^ ^ = Rf and vector spherical 
harmonics X^^^ = Y^^^ in terms of radial and angular coor- 
dinates r]± = r and ryy = {6>, ^}, respectively, and labeled by 
angular-momentum "quantum" numbers n = {^, m}. 

Because M^ and N^ form an orthonormal basis (due to the 
self-adjointness of the Helmholtz operator), the homogeneous 
photon DGFs G and C of Sec. II A can be expressed in such a 
basis as: 105-10^ 



G(/c; X, x') 



v±{^)v±{^') 



2ik 



(5fx - XM + V /^-'^^- W ^ ^n (XO + Xn,MN+(x) N" (x^ V±{^) > V±{^') 

^ ^ ^lXn,^M-(x)0M+(xO+Xn,MN-(x)0N+(xO v±{^)<v±{^V 



|: V X G, respectively, where the coefficients Xn are 

-) and 



andC 

determined by taking the Wronskian of the outgoing (- 

incoming (+) solutions. 

The SIE matrices appearing in Eq. 18 involve inner prod- 
ucts of Eq. 37 with basis functions f3n defined at the surfaces 
of the bodies. (Note that because the Green's functions are 
evaluated on the surface, inclusion of the delta- function term 



I 

is crucial^^^.) Just as the vector fields M^ and N^ form a 
convenient basis in which to expand waves propagating inside 
and outside dV, so too do the tangential components X^ and 
Zn form a suitable basis in which to express surface-current 
basis functions Pn defined on dV. In this case, as in the case 
of RWG basis functions,^^ (3^ can be chosen to be purely elec- 
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trie (E) or purely magnetie (M), so that 

Pn,E = ( Q^ I ' Pn,M = 



by: 







(38) 



Moreover, the orthogonaHty relations of the tangential vector 
fields, (X^,Xn) = (Z^,Z^) = dmn and (X^,Zn) = 0, 
mean that only basis functions with the same n and same po- 
larization have non-zero overlap. These surface currents form 
a complete basis and satisfy convenient orthogonality rela- 
tions with the corresponding vector fields: 



(X^,N±) = (Z^,M±)=0 



(39) 

(40) 
(41) 



with inner products (•, •) corresponding to surface integrals 
over the tangential coordinates evaluated at the surface dV, 
i.e. {(p^i/j) = ^^^d^ryy j7(7?x,'^7||)^*V^, where J denotes 
the Jacobian factor for the coordinate system. 

The combination of these orthogonality relations and the 
Green's function expression of Eq. 37, implies that the G 
matrices arising in the SIE formulation for interface dV will 
be diagonal and known analytically in this basis. Therefore, 
choosing this basis simplifies the calculation of <l>, as illus- 
trated in the next sections. 



2. Heat transfer formulation 



(G^^ G\ 
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where G±_ and G\\ = G±{E ^ JY) are 2 x 2 block matrices 



Gl 






.(42) 



Here, the subscripts ± and || refer to the two decoupled polar- 
ization states, corresponding to purely electric E and purely 
magnetic M surface currents, respectively. The separability 
of the two polarizations means that the flux spectrum ^ can 
be written in the form ^ = "^ ^p, with ^p denoting the 
contribution of the p polarization. From the definitions of the 
r functions, it follows that the two are related to one another 
by ^11 =^^{Z ^ 1/Z). 

In the subsequent sections, we derive semi-analytical ex- 
pressions for ^ in special geometries involving isolated and 
interacting plates and spheres. The symmetry of these geome- 
tries make it convenient to represent the SIE matrices using 
Fourier and spherical- wave surface basis functions, described 
in Appendix A. Our final expressions agree with previous for- 
mulae derived using the scattering-matrix approach.^'^"^'^^'^^ 



Expression of the homogeneous Green's function in the in- 
terior of each high- symmetry body r in the basis specialized 
for that body yields a block-diagonal self-interaction matrix 
G^ with matrix elements G^^^qq, = (/3;;,^g, T^ ^ /3; g,) - 
Smn, where Q denotes polarization. In contrast, the lack 
of any orthogonality relations between wave solutions con- 
structed for different, unrelated bodies means that the inter- 
action matrices G^'^^ are dense, i.e. the matrix elements 

g: 



0,1 
mn,QQ' 



(/3^ Q, F^ ^ P^^Qf) generally do not vanish. The 
outgoing fields into the intervening medium F^ ^ Pm,Q ^^^ 
to currents in body r are still known analytically from Eq. 37, 
described in terms of the wave solutions specialized to body r 
(albeit evaluated in the exterior medium 0), but in order to take 
inner products with /3^ g for a body r' we need to "translate" 
the solutions centered on r to the different basis of waves cen- 
tered on r' 7^ r. Such change of bases are often performed via 
"translation" and "conversion" matrices that are well-known 
and tabulated for most shapes of interest, ^^^ and immediately 
yield the SIE interaction matrices G^'^^ . 

For the remainder of the section, we restrict ourselves to sit- 
uations involving either a single body or two identical bodies 
described by a common set of basis functions, in which case 
the individual SIE matrices are block-diagonal in n and polar- 
ization. In particular, the G matrices for a given n are given 



3. Isolated plates 

We first consider the radiation of an isolated plate. Using 
the appropriate Fourier basis supplied in Appendix A 1 and 
the corresponding Green's function expansion of Eq. 37, the 
G± matrices for the plate are given by: 
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(43) 



where 7^ = y^l — (|k^|//c^)^ is the wavenumber in the z 
direction normalized by k^. It follows that the flux spectra for 
the two polarizations are given by: 



^j 



47r ' 
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Re 
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^11 



^^iz^ 



1 
Z^ 

(44) 



with Tr ^ = / 4;^ ^(k^) corresponding to integration over 
the parallel wavevector. Assuming a non-dissipative exter- 
nal medium (Im eq = Im /io = 0). and performing straight- 
forward algebraic manipulations, one obtains the well-known 
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formula for the emissivity of the plate: ^^ 

1 r ^'k^ 



^{uj) 



Sir J, {27ry 



Y, ep(k^,a;), (45) 



P={^,ll} 



where ep = ^(1 — |rpp) denotes the directional emissivity 
of the plate for the p polarization, expressed in terms of the 
Fresnel reflection coefficients:^^ 



70 71 

Zq Zi 

Zq Z\ 



r±{Z^ 
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(46) 



5. Two plates 

Finally, we consider the heat transfer between two paral- 
lel, semi-infinite plates separated by distance d. Just as in the 
case of isolated plates, it is convenient to express the G± ma- 
trices in the Fourier basis supplied in Appendix A 1 . Here, in 
addition to the self-interaction matrices 
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(52) 



for r = 
trices 



\ Zq / 

1,2, one obtains the interaction or "translation" ma- 



4. Isolated spheres 

We now consider the radiation of an isolated sphere. Using 
the appropriate vector spherical wave basis supplied in Ap- 
pendix A 2 and the corresponding Green's function expansion, 
the G± matrices for the sphere are given by: 
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where the exponential factors above couple or "translate" 
waves arising in different origins. Straightforward matrix al- 
gebra yields the following flux spectra for the two polariza- 
tions: 
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where f{z) = {l/z -\- d/dz)f, J£ and hi are Bessel func- 
tions of the first and second kind, respectively, and Zr = k^R. 
Employing a number of well-known properties of spherical 
Bessel functions, such as the Wronskian identity j[{z)hi{z) — 
h'^{z)ji{z) = i/ z'^, one arrives at the following flux spectra 
for the two polarizations: 
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with Tr $ = ^^ ^ ^£rn corresponding to a sum over the 
angular-momentum quantum numbers. Assuming vacuum as 
the external medium (sq = fiQ = 1) and a non-magnetic 
sphere (/ii = 1), one obtains the well-known formula for the 
emissivity of a sphere in vacuum:^^ 
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(51) 



where ni = y^ is the index of refraction of the sphere. 



where pp = |1— Tpfp 6^*^°^°^ p and r^ is the Fresnel reflection 
coefficient of plate q for the p polarization given in Eq. 46. As- 
suming a non-dissipative external medium (Im eq = Im //q = 
0), and performing straightforward algebraic manipulations, 
one obtains the well-known formula: ^^ 



^{Uj) = <3^prop(^) + ^evan(c^), 



with 
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{2^y 
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(58) 



where e^ denotes the emissivity of plate q for the p polar- 
ization, and where $ has been conveniently decomposed into 
far-field (propagating) and near-field (evanescently decaying) 
contributions. 



C. BEM discretization via RWG basis 

In contrast to spectral methods, BEMs discretize the sur- 
faces of the bodies into polygonal elements or "panels", and 
describe piecewise functions in each element by low-degree 
polynomials. ^^'^^ The most common BEM technique employs 
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0.1 
separation d/2R 

Figure 4. Ratio V, = H/aT^A of the heat-transfer rate H between 
gold spheres of radii R = 0.2/im and the Stefan-Boltzmann law 
aT^A, where A = AttR^ is the surface area of the spheres, with one 
sphere held at T = 300 K and the other held at zero temperature, as 
a function of their surface-surface separation. (Inset:) Flux spectra 
^(u) per unit area A (units of /xm^), of the two spheres at d = i? 
(green circles) and of an isolated sphere (blue circles). 



a so-called RWG basis of vector- valued polynomial functions 
defined on a mesh of triangular panels. ^^ Such a basis is ap- 
plicable to arbitrary geometries and yields results that con- 
verge with increasing resolution (smaller triangles), where 
variants with different convergence rates depend upon the de- 
gree of the polynomials used in the triangles (which can be 
curved). The simplest discretizations involve degree- 1 poly- 
nomials and flat triangles, where the error decreases at least 
linearly with 1 /diameter of the triangles, but can converge 
faster with adaptive mesh refinements.^^ In contrast to spectral 
methods, the Gnm integrals here must be performed numeri- 
cally and the resulting G matrices are dense, but thankfully 
fast techniques to perform these integrals are well established 
and need only be implemented once for a given RWG basis, 
independent of the geometry. ^"^'^^'^-^ One such implementation 
is the free-software solver SCUFF-em,^^^ which we exploit in 
this section to compare results from BEM RWG to known re- 
sults for spheres; the same code is employed in Sec. V to ob- 
tain results in new and more complex geometries. 

The heat-transfer rate H between two spheres was recently 
obtained numerically by Ref. 29. In contrast to scattering- 
matrix methods or the FSC formalism above, the method 
of Ref. 29 involves straightforward integration of the inho- 
mogeneous Green's function of the geometry over the vol- 
umes of the two spheres, expressed in terms of a special- 
ized spherical-wave basis expansion with coefficients deter- 
mined by enforcing continuity of the fields across the vari- 
ous interfaces. The result of the integration is an exponen- 
tially convergent semi-analytical formula of the kind derived 
in Sec. IV B. Figure 4 compares the results of the BEM 



RWG method (red circles) against those obtained by evalu- 
ating the semi-analytical formula of Ref. 29, truncated at a 
sufficiently large but finite order (solid lines). In particular, 
the heat transfer ratio H = H/crT^A is plotted as a func- 
tion of surface-surface separation d for gold spheres of radius 
R = 1/im, where one sphere is held at T = 300 K while 
the other is held at zero temperature, and where crT^A is the 
Stefan-Boltzmann (SB) law (for a planar black body), with 
(T = tt'^ k% / {60h^ c^) and A the surface area of the spheres. 
The inset of the figure also shows the corresponding flux spec- 
tra ^ of both interacting (d = R) and isolated spheres, nor- 
malized by A and plotted over relevant wavelengths A > At, 
where At = hc/ksT ^ 7.6/im denotes the thermal wave- 
length corresponding to the peak of the thermal spectrum. In 
both cases the BEM results (circles) are shown to agree with 
the corresponding semi-analytical formulas (in the case of iso- 
lated spheres, the flux spectrum is compared against Eq. 51). 



V. APPLICATIONS 

In this section, we illustrate the generality and broad appli- 
cability of the FSC formulation by applying the BEM RWG 
method to obtain new results in complex geometries. As dis- 
cussed above, most calculations of heat transfer have focused 
primarily on semi-infinite planar bodies. ^^ Finite bodies only 
recently became accessible with the development of sophis- 
ticated spectral methods, ^^~^^'^^ albeit for highly symmetric 
bodies with smooth shapes (e.g. spheres) for which conve- 
nient spectral bases exist. Here we will focus instead on ge- 
ometries involving finite bodies with sharp corners (combina- 
tions of finite plates, cylinders, and cones) that pose no chal- 
lenge for the BEM RWG method but which prove difficult to 
model via spectral methods. 



A. Plates and cylinders 

To begin with, we extend the calculation of heat transfer 
between planar semi-infinite plates to the case of finite plates, 
which quantifies the influence of lateral size effects in that 
geometry. Figure 5 shows the ratio H = H/crT^A of the 
heat-transfer rate H between finite circular plates of thick- 
ness L and lateral size 2R and the SB law, with one plate 
held at T = 300 K and the other held at zero temperature, 
as a function of their surface-surface separation d. H is plot- 
ted for multiple aspect ratios 2R/L (solid circles), with fixed 
L = 0.2/im. For comparison, we also plot the heat- transfer 
ratio Hoc for semi-infinite (R -^ oo) plates of the same thick- 
ness (black solid line), which is obtained analytically from 
the absorptivity of the plates via Kirchhoff 's law of thermal 
radiation.^'^^^ As expected, one finds that at small d, near-field 
effects dominate and H ~ 1/V^ for both finite and semi- 
infinite plates. In contrast, at asymptotically large d, the fi- 
nite plates behave like dipoles and one finds that H ~ l/d^, 
whereas the semi-infinite transfer rate approaches a constant 
l-Loo{d -^ oo) <C 1 independent of d; the rate is signifi- 
cantly smaller than that of a perfect black body because gold is 
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separation d/L 

Figure 5. Ratio H, = H/aT^A of the heat-transfer rate H from a 
finite, gold circular plate of lateral size 2R and thickness L = 0.2/im 
held at T = 300 K, to an identical plate held at zero temperature, 
and the SB law aT^A (where A = ttR^ is the "interaction" surface 
area of the plates), as a function of their surface-surface separation 
d. 7i is plotted for multiple aspect ratios 2R/L (circles). The solid 
black line corresponds to the heat-transfer ratio Hoo = Hoo{R -^ 
oo) obtained upon taking the limit R ^ oo, which is computed via 
the semi-analytical formula in Ref. 110. (Inset:) Heat-transfer rate 
between two plates at a fixed separation d = O.IL (solid circles), 
and heat radiation of an isolated plate (open circles) or sphere (thick 
solid line), as a function of lateral size or diameter. 
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Figure 6. Flux spectrum ^{uj) of an isolated gold cylinder of length 
L and radius R = 0.2/im, normalized by its corresponding surface 
area A, and plotted for multiple aspect ratios L/R (solid circles). 
The solid line shows $ in the limit L ^ oo of a semi-infinite cylin- 
der, as computed by the semi-analytical formula of Ref. 34. (Inset:) 
Heat-transfer ratio 7i of heat transfer H from a room-temperature 
cylinder of aspect ratio L/R = 5 to an identical cylinder at zero 
temperature, and the SB law aT^A, with A = 2tiR{R + L) de- 
noting the total surface area of each cylinder, as a function of their 
surface-surface separation d. H is plotted for both parallel (0 = 0) 
and crossed (0 = 90°) cylinder configurations, with the shaded re- 
gion corresponding to intermediate 0. 



highly reflective. As i? ^ oo, the BEM results approach Hoo 
for all separations d, albeit at different rates, where smaller 
separations converge faster than larger separations. 

To quantify finite-size effects, the inset of Fig. 5 shows 
V,/l-Loo for isolated and interacting plates (at a single sepa- 
ration d = O.IL) as a function of R. As above, in the limit of 
large R :^ Xt ^ L, such that the dominant wavelengths and 
corresponding skin depths S = c/lm. ^Jsijj are much smaller 
than the lateral dimensions of the plates, H -^ Hoc- In the 
case of isolated plates, the relevant lengthscales are At and ^, 
whereas in the case of interacting plates, the separation d also 
factors into the convergence rate: the increasing contribution 
of (long-wavelength) near-field effects to the heat transfer at 
smaller separations means that smaller separations converge 
faster to the Hoc result than larger separations. (For the par- 
ticular separation (i = O.IL plotted here, near field effects are 
large enough to cause the convergence rate of the interacting 
plates to be significantly larger than that of the isolated plate.) 
At intermediate i? < L <C At, the plates no longer resem- 
ble plates but rather elongated cylinders, leading to significant 
deviations in H. 

Compared to the heat radiation of semi-infinite cylinders 
{LjR — > oo for fixed K), studied previously by Ref. 34, 
the radiation of finite cylinders displays a number of inter- 
esting features. (Note that H here includes radiation emitted 



in both the axis-parallel, Hy, and -perpendicular, H^, direc- 
tions.) First, due to the finite value of L, in the limit R ^ oo 
the radiation of our finite cylinders is best characterized by 
the radiation of thin plates with Hy ^ 1-L_\_. Not surprisingly, 
we find that H ^ Hoc from below sls R -^ oo, in contrast 
to what is observed in the semi-infinite case where Hoo is 
approached from above. ^^ Second and most interestingly, we 
find that below a critical R, determined by the smallest skin- 
depth S ^ 20nm of Au over the relevant thermal wavelengths, 
the radiation normalized by surface area increases with de- 
creasing R, leading to non-monotonicity. Such behavior is 
unusual in that in this R ^ S regime, bodies most often be- 
have like volume emitters, causing H to grow with the vol- 
umes rather than surfaces of the bodies (as observed in the 
case of semi-infinite cylinders). ^^ Indeed, we find that for di- 
electric bodies with small and positive £, one obtains the usual 
volume-dependence of H. In contrast, the enhancement in 
Fig. 5 arises because for small R, the cylinders act as metallic 
dipole emitters, whose radiation is increasingly dominated by 
H|| as i? ^ and whose quasistatic (long wavelength) paral- 
lel polarizability grows with decreasing R (a consequence of 
the increasing anisotropy of the cylinder and large Im e)}^^'^^^ 
For sufficiently small R, the heat transfer per unit area of the 
uniaxial cylinders can greatly exceed that of the semi-infinite 
plate, i.e. H ^ Hoo- (The dipole model also predicts that H 
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will eventually vanish sls R -^ 0, but only at radii too small 
to be easily calculated by BEM. We intend to explore these 
phenomena more fully in subsequent work.) 

It is also interesting to study the convergence of the cylin- 
der radiation rate with L, comparing our results against the 
semi-analytical results obtained in the special case of semi- 
infinite (L ^ oo) cylinders. ^"^ We also consider the heat trans- 
fer between non-uniaxial (parallel) cylinders. Figure 6 shows 
the flux spectra <l> of isolated cylinders of radius R = 0.2/im 
and varying lengths L; for comparison, we also plot the spec- 
trum of the semi-infinite cylinders^"^ (solid lines). As before, 
<l> is normalized by the surface area A of each object. (For 
the relevant wavelength range shown in the figure, R is sev- 
eral times S, which means that most of the radiation is com- 
ing from sources near the surface of the objects. ^'^) We find 
that for L/R ^ 2 (not shown), corresponding to nearly- 
isotropic cylinders, $ is only slightly larger than that of an 
isolated sphere due to the small but non-negligible contri- 
bution of volume fluctuations to <l>. As L/R increases, ^ 
increases over all A, and converges towards the L ^ oo 
limit (black solid line) as A -^ 0, albeit slowly. Moreover, 
^L ^ ^oo at particular wavelengths, a consequence of geo- 
metrical resonances that are absent in the semi-infinite case — 
away from these resonances, <l> clearly straddles the L ^ oo 
result so long as A < L. As in the case of finite plates, 
the ^ of interacting cylinders exhibits significant enhance- 
ment at large A due to near-field effects, so that 1-L ^ oo 
with decreasing separation d. The enhancement is evident 
in Fig. 6, which shows H over a wide range of d for both 
parallel- {9 = 0) and crossed-cylinder (0 = 90°) configura- 
tions, with one cylinder held at T = 300 K and the other at 
zero temperature (both cylinders have aspect ratio L/R = 5). 
We find once again that there are two very distinct separation 
regimes of heat transfer: at large d^ R, the cylinders act like 
dipole emitters and 1-L/l-Loo ^ ^/d^ <^ 1 whereas at small 
d <^ R, flux contributions from evanescent waves dominate 
and 1-L/l-Loo ^ ^/\fd ^ 1. Comparing the heat transfer H in 
the parallel and crossed-cylinder configurations, we find that 
H\\/H±_ ^ 1 at large d ^ R but increases significantly at 
smaller d <^ R, again due to near-field effects: in the d ^ 
limit, H is dominated by closest surface-surface interactions, 
so H\\/Hj_ ~ L/R -^ 5. As expected, H\\/H± ^ oo as 
L ^ oo because the increased "interaction" area in this limit 
favors the parallel over the crossed configuration. Specifically, 
whereas H grows linearly with L in the parallel configuration, 
it grows sub-linearly (and asymptotes to a finite value in the 
L ^ oo limit) in the crossed configuration due to the dimin- 
ishing contributions of near-field and radiative flux between 
surface elements in the extremities of the cylinders. 



B. Cones 

Finally, motivated by recent predictions,^^ we consider the 
heat transfer between finite cones and plates. In Ref. 69, 
the cone-plate geometry (with a semi-infinite plate) was ob- 
tained using a "hybrid" scattering-BEM method^^ based on 
the scattering-theory formulation of Ref. 32. (In contrast to 
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Figure 7. Heat- transfer rate H from a room- temperature gold cone of 
base radius R = l/xm and length L = 2R, to either an identical cone 
(green circles) or a gold plate of radius R and thickness h = 0.2/im, 
held at zero temperature, as a function of their surface-surface sep- 
aration d. H is normalized by the Stefan-Boltzmann law aT^A, 
where A is the surface area of the cone. (Inset:) Flux spectrum ^(cj) 
of the cone-plate configuration at a single separation d = 0.2L, nor- 
malized by the area of the cone. The two surface-contour plots show 
the distribution of flux pattern on the surfaces of the bodies at two 
wavelengths, A ^ SOL and A ^ 2.2L, where white/black denotes 
the maximum/minimum flux at the corresponding wavelength. 



semi-infinite plates or spheres, the scattering-matrix of a cone 
cannot be easily obtained analytically, and was instead com- 
puted numerically by exploiting the BEM method in combi- 
nation with a multipole basis of cylindrical waves.) Here, 
in addition to extending these predictions to the case of fi- 
nite plates, we consider the heat- transfer rate between two 
oppositely oriented cones. Figure 7 shows the heat- transfer 
rate h (as in the previous section, h = H/cjT^A where here 
A = ttR^ is the projected area of the cone) from a cone of 
radius R = 0.5/im and length L = 2R to either an iden- 
tical cone rotated by 180° (green circles) or a plate of ra- 
dius R and thickness L = 0.2/im (red circles), as a func- 
tion of their surface-surface separation d. As before, we con- 
sider gold bodies, with one held at 300 K while the other is 
held at zero temperature Similar to Ref. 69, we find that the 
heat transfer rate H ~ log{d) varies logarithmically with d 
at short separations d <C L <C At, a consequence of near- 
field interactions and the finite size of the cone.^^"^ While h 
exhibits similar scaling with d for both geometries, h turns 
out to be almost two orders of magnitude smaller at small 
d <C L in the cone-cone geometry, as would be expected 
from a proximity-approximation (PA) model. ^^^ The situation 
is reversed at large separations d^ \t ^ L\ beyond a criti- 
cal d ^ 7L, the cone-cone heat transfer becomes larger than 
the cone-plate transfer. The reversal is expected on the basis 
that at these separations, the two bodies act like fluctuating 
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dipoles oriented mainly along their largest dimension (along 
the axis of symmetry for the cone and along the lateral dimen- 
sion for the plate), in which case the cone-plate interaction re- 
sembles the interaction of two orthogonal dipoles whereas the 
cone-cone interaction resembles the interaction of two par- 
allel dipoles. Another interesting feature of the heat-transfer 
in this geometry is that the spatial distribution of pattern over 
the plate exhibits a local minimum directly below the tip of the 
cone, a consequence of the dipolar field induced on the cone at 
long wavelengths.^^ Here, we observe a similar phenomenon 
but we find that the finite size of the plate significantly alters 
the scope of the anomalous radiation pattern. In particular, 
whereas Ref . 69 found this effect to persist over a wide range 
of wavelengths (surviving even in the total or integrated radi- 
ation pattern), we find that in the finite-plate case it disappears 
much more rapidly with decreasing wavelength. 



VI. CONCLUSION 

The FSC approach to non-equilibrium fluctuations pre- 
sented here permits the study of heat transfer between bod- 
ies of arbitrary shape, paving the way for future exploration 
of heat exchange in micro structured geometries that until now 
remain largely unexplored in this context. Our formulation 
shares many properties with previous scattering-matrix for- 
mulations of radiative heat transfer, e.g. our final expressions 
involve traces of matrices describing scattering unknowns, but 
differs in that our "scattering unknowns" are surface currents 
defined on the surfaces of the bodies rather than incident and 
outgoing waves propagating into and out of the bodies. ^^~^^ 
As argued above, this choice of description has important 
conceptual and numerical implications: it allows direct ap- 
plication of the surface-integral equation formalism as well 
as the boundary-element method. When specialized to handle 
high- symmetry geometries using special functions that exploit 
those symmetries, our approach can be used to obtain fast- 
converging semi-analytical formulas in the spirit of previous 
work based on spectral methods. ^^'^"^'^^ Moreover, it can also 
be applied as a brute-force method, taking advantage of ex- 
isting, well-studied, and sophisticated BEM codes (with no 
modifications), to obtain results in arbitrary/complex geome- 
tries. 

While the main focus of this work was on exploring some 
of the ways in which the FSC formulation can be applied 
to study non-equilibrium heat transfer, we believe that anal- 
ogous techniques can be used to derive corresponding FSC 
approaches to other fluctuation phenomena, including near- 
field fluorescence, ^^^ quantum noise in lasers, ^^^ and non- 
equilibrium Casimir forces,^^'^^^ an idea we plan to explore 
in future work. Furthermore, although our calculations here 
focused on geometries involving compact bodies, the same 
heat-transfer formulas derived above apply to geometries in- 
volving infinitely extended/periodic bodies (of importance in 
applications of heat transfer to thermophotovoltaics). Modi- 
fying BEM solvers to handle periodic structures, however, is 
non-trivial^ ^^~^^^ and we therefore consider that case in a sub- 
sequent publication. 



Finally, although Eq. 18 is already well-suited for efficient 
numerical implementation, its computational efficiency may 
be improved by adopting a modified formulation in which the 
dense G matrices are replaced by certain sparse matrices in- 
volving overlap integrals among basis functions. In addition 
to reducing the computational cost of the trace in Eq. 18, this 
approach has the advantage of allowing the computation of 
other fluctuation-induced quantities such as non-equilibrium 
Casimir forces and torques. This alternative formulation will 
be discussed in a forthcoming publication. 
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Appendix A: Eigenf unctions of the Helmholtz equation 

In this section, we provide and exploit the standard Fourier 
and spherical-wave eigenfunctions of the vector Helmholtz 
operator, obtained by solving Eq. 34 in planar and spherical 
coordinates, ^^^'^^^ to obtain the coefficients x ^nd n appearing 
in the Green's function expansion and orthogonality relations 
of Eq. 37 and Eqs. 39-41, respectively. 



1. Fourier basis 

In planar geometries, described by normal and tangential 
coordinates z and x^ , respectively, the eigenfunctions of the 
Helmholtz operator, labeled by Fourier wavevectors k^ per- 
pendicular to the z axis, are given by: 



N^fc (^,x^)=^±(fc,z) 



TyZk^(Xi) + -^e z 



where (p^^^^ = ^^e^'^^-^±±ik,^ ^ k^ = ^Jk^ - |k_L|2, and 
where the tangential fields Xkj^ and Zkj^ = z x Xkj^ are: 



Xkjx^) = 



ki 



(zxki)e*''-^ 



ki x I 



The precise form of the Fourier functions (j)^ 
pends on whether one desires a solution involving outgoing 
(+) or incoming (— ) fields, or equivalently, fields propagat- 
ing away or toward the origin. The corresponding x ^^^ ^ 



(Al) 

(A2) 
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coefficients appearing in the Green's function expansion and 
orthogonaHty relations are given by: 



<.fc..sw = <t>^ik^^) 


(A3) 


4^,fe.,MW = T7'^*(fc^2) 


(A4) 


Xk,,fe. - 2^^ 


(A5) 



with -f = k^/k = v^l- |k^|7(6:/icj2). 



2. Spherical multipole basis 

In spherically symmetric geometries, described by nor- 
mal and tangential coordinates r and {6>, 0}, respectively, the 
eigenfunctions of the Helmholtz operator, labeled by angular- 
momentum quantum numbers i and m, are given by: 



Mf^{r,0,^) = Rf{kr)X,m{OA)^ 



Nf^{r,0,^) = Rf{kr)Z, 



i{i^l) 



Rf{kr)Yim{0,^)v, 



?± 



where R^ and Yim denote spherical Hankel functions and 
spherical harmonics, ^^ respectively, and where the tangential 



fields X^^ 
are: 



V^(^+i; 



:(r X Vjy^^andZ, 
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Above, we defined f{z) = {l/z -\- d/dz)f for brevity. The 
precise form of the spherical Bessel radial function 



Rf — 



.(1) 



\Ji 



depends on whether one desires a solution corresponding to 
outgoing (+) or incoming (— ) waves toward the origin, or 
equivalently, a solution that is well-behaved at the origin or 
at infinity. The x and tz coefficients appearing in the Green's 
function expansions and orthogonality relations are given by: 



4.m.Ei'r)='r^Rfikr) 



'«tm,MW=«^^^^*(^0 



Xlr. 
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Appendix B: Equivalence principle 

In this section, we provide a compact derivation and re- 
view of the equivalence principle of classical electromag- 
netism (closely related to Huy gens' principle^^), which ex- 
presses scattered waves in terms of fictitious equivalent cur- 
rents in a homogeneous medium replacing the scatterer.^"^ The 



equivalence principle is usually derived in a somewhat cum- 
bersome way from a Green 's-function approach,^"^'^^ but a 
much shorter proof can be derived from the differential form 
of Maxwell's equation. Understanding this result is central to 
our FSC formulation of heat transfer. 

As before, we restrict ourselves to linear media, for which 
Maxwell's equations can be written as: 



d_ 
dt 



[0 + X ^ 0] 



(Bl) 



M 



with x^ denoting convolution with a 6 x 6 susceptibility tensor 



/i 



Consider an arbitrary incident wave (j) which solves the 
source -free Maxwell's equations in some x medium with no 
current sources: Mcj) = ^[0 + X^^]- The equivalence princi- 
ple states that, given any arbitrary hut finite domain V, one can 
always choose an equivalent surface current ^ that generates 
the same incident field (j) inV. To show that such a surface 
current exists, define the field 



elsewhere 



(B2) 



It follows that ^ satisfies the source-free Maxwell's equations 
in both the interior and exterior regions — the only question is 
what happens at the interface dV . In particular, the disconti- 
nuity of (}) at dV produces a surface delta function Sqy in the 
spatial derivative M0, and so in order to satisfy Maxwell's 
equations with this ^ one must introduce a matching delta 
function, a surface current ^, on the right-hand side. [Here, 
5dv is the distribution such that ^(5ay0(x) = ^^^ 0(x) 
for any continuous test function 0.] Specifically, letting n be 
the unit inward-normal vector ^^^, only the normal derivative 
n • V contains a delta function (whose amplitude is the mag- 
nitude of the discontinuity), which implies a surface current: 



e = (e(/>)fey = 



nxH 

-n X E 



^d\ 



(B3) 



where B is the real- symmetric unitary 6x6 matrix: 



6 



nx 



6- 



6^ 



e*. 



(B4) 



That is, there is a surface electric current given by the surface- 
tangential components n x H of the incident magnetic field, 
and a surface magnetic current given by the components 
— n X E of the incident electric field. These are the equivalent 
currents of the principle of equivalence (derived traditionally 
from a Green's function approach^"^'^^ and from which Huy- 
gens's principle is derived^^). 
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1. Application to surface integral equations 

The equivalence principle is of fundamental importance to 
SIE formulations of EM scattering. Consider two regions 
and 1, described by volumes V^ and V'^ and susceptibilities 
X^ and x^^ respectively, separated by an interface dV^. As 
before, one can express the total fields (j)^ = ^^++^^~ in each 
region r, in terms of incident 0^+ and scattered 0^~ fields. 
The principle of equivalence describes an equivalent, fictitious 
problem, involving fields 

elsewhere 



need only prove reciprocity and definiteness of the homoge- 
neous Green's function (trivial to show in that case since the 
homogeneous DGF is known analytically), here we consider 
the more general case of inhomogeneous media. Reciprocity 
is well known,^^'^^^"^^^ and positivity follows from general 
physical principles (currents always do nonnegative work in 
passive materials^^'^"^'^^^'^^^), but our goal here is to derive 
them using the same language employed in our derivations 
above. More specifically, we explain the source of the sign- 
flip matrices S and S, which often go unmentioned because 
many authors consider only 3x3 Green's functions (relating 
currents to fields of the same type). 



and surface currents ^ = 60^ = 6^^ at the dV^ interface, 
where the second equality follows from continuity of the tan- 
gential fields. Since = in F^, it follows that one can re- 
place x^ with any other local medium and yet ^ still satisfies 
Maxwell's equations. In particular, replacing x^ with x^ im- 
plies that one can write the scattered field (j)^~ =T^ ik^inV^ 
as the field produced by the same fictitious surface currents ^ 
in an infinite medium 0, with F^ denoting the homogeneous- 
medium Green's function of the infinite medium. 

A similar argument applies if one is interested in the field 
in medium 1, except that the sign of the fictitious currents 
is reversed to — ^ in order to account for the direction of the 
discontinuity in going from 1 to in this case. In particular, 
one can write the scattered field (j)^~ = —T^i^^inV^ as the 
field produced by a fictitious surface current — ^ in an infinite 
medium 1. 



Appendix C: Reciprocity and definiteness 

In this section, we present a brief review of the reciprocity 
relations and definiteness (positivity) properties of the DGF, 
F, connecting surface currents ^ to fields (j) = F ^ ^, in dis- 
sipative media, and explain how these relate to corresponding 
properties of the SIE matrices above (crucial to our deriva- 
tion of heat transfer in Sec. II). Although for our purposes we 



1. Green's functions 

It is actually easier to derive the reciprocity and definite- 
ness properties of F from the properties of L = (F^)~^, 
the Maxwell operator that connects fields (j) to currents ^ = 
L(j), because L is a partial-differential operator that can be 
written down explicitly starting from the (frequency-domain) 
Maxwell equations V x E = iufiH — M, V x H = 
—iuje'E + J, in terms of the permittivity £:(x, uj) and perme- 
ability ii{^^uj) tensors and electric J and magnetic M cur- 
rents. Specifically, the Maxwell operator 



L 



iuje Vx 
— Vx iujfi 



(CI) 



is neither complex-symmetric, Hermitian, anti- symmetric, 
nor anti-Hermitian in general. Using our previous definition 
of the inner product: 



= /e* E' + H* H^ 



E' 



it follows that the 6>j^-diagonal part of L is anti-Hermitian: 



Vx 



-Vx 



E' 



E* . V X H' - H* . V X E' 



/(V X E*) • H' - (V X H)* • E' = /- ( 



Vx 



-Vx 



E' 



where we have used the self-adjointness of Vx and assumed 
boundary conditions such that the ^ E* x H' + E' x H* 
boundary terms at infinity (from the integration by parts) van- 
ish. This is commonly attained by assuming loss in the materi- 
als so that the fields decay exponentially at infinity (assuming 
localized sources), or by imposing outgoing-radiation bound- 
ary conditions on F ^ at infinity. ^^ 



r 



Instead, reciprocity relations are normally derived for the 
unconjugated inner product: 



= /e'^-e' + h^-h', 



E' 
H' 



(C2) 
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under which the off-diagonal terms in L are still anti- 
symmetric while the diagonal terms are complex- symmetric, 
assuming reciprocal materials: e^ = e and 11^=11 (usu- 
ally the case except for magneto-optical and other more ex- 
otic materials^^'^^'^^^). Here, the transpose L^ of the oper- 
ator L means the adjoint of L under the unconjugated inner 
product (0, L(j)') = (L^0, (j)'). In order to make L fully sym- 
metric, it suffices to flip the sign of the magnetic components 
H ^ — H, an operation that can be expressed as a (real, self- 
adjoint, unitary) sign-flip matrix: 



S 



= s- 



:5' =5* 



(C3) 



That is, LS is complex- symmetric: (LS)^ = SL^ = LS, or 
equivalently, L^ = SLS = SLS~^. It follows that, 

{Ti.f = {L-^f = (L^)-^ = S{Ti.)S (C4) 

Alternatively, 

{^, r*cP') = 1 1 d3xdV'^^(x)r(x, y)(/.'(y), (C5) 

SO by inspection (r(x, y)^)^ = r(y,x)^^: transposing F^ 
corresponds to interchanging sources and fields. Thus, we ob- 
tain the reciprocity relation: 



r^^ = 5(r^)5, 



(C6) 



i.e. one can interchange sources and fields if one flips the sign 
of both magnetic currents and magnetic fields. 

We also expect the operators L and F^ to be negative- 
semidefinite on physical grounds, since —^{(j)^L(j)) = 

— I ((/), ^) = — ^ (F ^ ^, ^) is exactly the time-average power 

— |/E*-J + H*-M expended by the currents, which must be 
> in passive materials. ^"^ Indeed, one can show this directly. 



since the anti-Hermitian property of the off-diagonal part of L 
means that 



sym L = uj 



- Ivue 



- Im/i 



for isotropic materials. But both ujlvae and cj Im/i are > 
for passive materials (no gain).^^'^"^ Thus, it follows that L is 
negative-semidefinite, and so is L~^ = F ^ . 



2. SIE matrices 

The SIE matrices M = W~'^ are formed from a sum M. 
of Green's function operators F^^ in homogeneous regions 
r, expanded in a (real vector- valued) basis Pn by a Galerkin 
method, so that Mmn = (/^m, Mpn) = (/3m, Mpn). For any 
Galerkin method, it is easy to show that if M is self- adjoint 
or complex-symmetric, then M has the same properties. Sim- 
ilarly, any definiteness of A^ carries over to M. From the pre- 
vious section, since F^ is negative-semidefinite in any passive 
medium, it follows that any sum M of F^^ is also negative- 
semidefinite, and hence M is negative-semidefinite (sym M 
is Hermitian negative-semidefinite). 

As above, reciprocity requires some sign flips: M^ = 

SMS, so {M^)mn = Mnm = {Pn^M^m) = 

iPm.M^Pn) = iPm.SMSpn) = {S (3m. MS (3^). Fur- 
thermore, suppose that we use separate basis functions (3^ 
for magnetic currents and f3^ for electric currents, as is typ- 
ically the case in BEM (e.g. for an RWG basis^^'^^), so 
that Sf3^ = +/3^ and Sf3^ = -f3^. That is, we write 
currents as ^ = T^^Xn(3n = Y^x^p^ + ^^ Pn ^ so that 
5*^ = 'Yx^/3^ — x^ /3^ corresponds to a linear transfor- 
mation S on X that flips the sign of the x^ components. It 
follows that 



M^ = SMS. 
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